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Abstract: The purpose of this work is simulation of magnetised plasmas in the 
ITER project framework. In this context, Vlasov-Poisson like models are used 
to simulate core turbulence in the tokamak in a toroidal geometry. This leads 
to heavy simulation because a 6D dimensional problem has to be solved, 3D 
in space and 3D in velocity. The model is reduced to a 5D gyrokinetic model, 
taking advantage of the particular motion of particles due to the presence of a 
strong magnetic field. However, accurate schemes, parallel algorithms need to be 
designed to bear these simulations. This paper describes a Hermite formulation 
of the conservative PSM scheme which is very generic and allows to implement 
different semi-Lagrangian schemes. We also test and propose numerical limiters 
which should improve the robustness of the simulations by diminishing spurious 
oscillations. We only consider here the 4D drift-kinetic model which is the back- 
bone of the 5D gyrokinetic models and relevant to build a robust and accurate 
numerical method. 
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Test of some numerical limiters for the 
conservative PSM scheme for 4D Drift-Kinetic 

simulations. 

Resume : Ce travail concerne la simulation de plasmas magnetises dans le 
cadre du projet ITER. Pour cette application, des modcles de type Vlasov- 
Poisson sont utilises pour simuler la turbulence a coeur dans un tokamak, en 
geometrie toroidale. Ces etudes menent a resoudre des problemes dans un espace 
a 6 dimensions, 3D cn espace 3D en vitessc, qui sont tres lourds a simuler en 
terme de ressources informatiques. Le modele est reduit a un modele gyrocine- 
tique 5D en exploitant les caracteristiques de ce plasma, dont le mouvement 
des particules est fortement influence par la presence d'un champ magnetique 
intense. Cependant, il est necessaire de mettre au point des schemas precis et 
des algorithmes paralleles pour mener ces simulations. Ce rapport decrit une 
formulation de type Hcrmitc du schema conservatif PSAT qui est tres generiquc 
et qui permet d'implementer different schemas semi-Lagrangiens. Nous testons 
et proposons egalement des limiteurs numeriques de pente qui doivent accroitre 
la robustessc des simulations cn rcduisant les oscillations d'origine numerique. 
Dans ce travail, nous I'utilisons pour resoudre le modele drift-kinetic 4D, qui 
est le squelette du modele gyrocinetique 5D. Ce modele 4D est sufiisamment 
pertinent pour la conception d'une methode numerique robuste et precise pour 
le modele 5D 

Mots-cles : simulation numerique, schema conservatif, ITER, turbulence 
plasma 
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1 Introduction 

The ITER device is a tokamak designed to study controlled thermonuclear fu- 
sion. Roughly speaking, it is a toroidal vessel containing a magnetized plasma 
where fusion reactions occur. The plasma is kept out of the vessel walls by a 
magnetic field which lines have a specific helicoidal geometry. However, turbu- 
lence develops in the plasma and leads to thermal transport which decreases the 
confinement efficiency and thus needs a careful study. Plasma is constituted of 
ions and electrons, which motion is induced by the magnetic field. The charac- 
teristic mean free path is high, even compared with the vessel size, therefore a 
kinetic description of particles is required, see Dimits [5]. Then a full 6D Vlasov- 
Poisson model should be used for both ions and electrons to properly describe 
the plasma evolution. However, the plasma flow in presence of a strong mag- 
netic field has characteristics that allow some physical assumptions to reduce 
the model. First, the Larmor radius, i.e. the radius of the cyclotronic motion of 
particles around magnetic field lines, can be considered as small compared with 
the tokamak size and the gyration frequency very fast compared to the plasma 
frequency. Thus this motion can be averaged (gyro-average) becoming the so- 
called guiding center motion. As a consequence, 6D Vlasov-Poisson model is 
reduced to a 5D gyrokinetic model by averaging equations in such a way the 6D 
toroidal coordinate system (r, 0, 0, Ur, we, w^) becomes a 5D coordinate system 
(r, ^, (/), Wjj , /i), with W|| the parallel to the field lines component of the velocity 
and ^ = m v\/2B the adiabatic invariant which depends on the norm of the 
perpendicular to the field lines components of the velocity v]_ , on the magnetic 
field magnitude B and on the particles mass m. Moreover, the magnetic field is 
assumed to be steady and the mass of electrons rUe is very small compared to 
the mass of ions mj. Thus the cyclotron frequency w^.e — Qi^e B/rrii^e is assumed 
to be much higher for electrons than for ions > > Wi. Therefore the electrons 
are assumed to be at Boltzmann equilibrium, i.e. the effect of the electrons 
cyclotronic motion is neglected. The 5D gyrokinetic model then reduces to a 
Vlasov like equation for ions guiding center motion; 

where f^{X,v\\) is the ion distribution function with X — {r,6,(l)), velocities 
dX/dt and dv\\/dt define the guiding center trajectories. 
If V(x,t;||) • {dX/dt,dv\\/dtY = 0, then the model is termed as conservative. 
This equation for ions is coupled with a quasi-neutrality equation for the electric 
potential ^(i?) on real particles position, with R = X — (with the Larmor 
radius) : 

--^Vi • (noVi$) + ■^(*- < $ >^) = / Ud^ldv\\ - no (2) 

where no is the equilibrium electronic density, Tg the electronic temperature, 
e the electronic charge, k the Boltzmann constant for electrons and oji the cy- 
clotronic frequency for ions. 

These equations are of a simple form, but they have to be solved very ef- 
ficiently because of the 5D space and the large characteristic time scales con- 
sidered. However, the adiabatic invariant p. acts as a parameter, thus it could 
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easily be parallelized. Moreover, we can see that for each /z, we have to solve 
a 4D advection equation, as accurately as possible but also taking special care 
on mass and energy conservation, especially in this context of large character- 
istic time scales. The maximum principle that exists at the continuous level 
for the Vlasov equation should also be carefully studied at discrete level be- 
cause there is no physical dissipation process in this model that might dissipate 
over/undershoots of the scheme. These studies will be achieved first on a rele- 
vant reduced model, the 4D drift-kinetic model which corresponds to ([T]) with 
fj, = instead a range of /i values (theoretically R+).This work follows those 
of Grandgirard et al in the GYSELA code, see [S] and [S]. The geometrical 
assumptions of this model for ion plasma turbulence are a cylindrical geometry 
with coordinates (r, 9, z, v\\) and a constant magnetic field B = Cz, where Cz 
is the unit vector in z direction. In this coUisionless plasma, the trajectories are 
governed by the guiding center (GC) trajectories: 

dr dO dz dv\\ qi 

with VGC = {E X B)/B'^ and E = -V$ with $ the electric potential. 

The Vlasov equation governing this system, with the ion distribution function 

/(r, 0, z, W|| , i), is the following: 

dtf + VGc.drf + VGCedef + vi\dj + —Ezd^J = 0. (4) 

This equation is coupled with a quasi-neutrality equation for the electric po- 
tential $(r, 0, z) that reads the same as for the 5D gyrokinetic model ^ with 
fi = 0. 

Let us notice that the 4D velocity field a = (ugc"^, VfjCj, , tiy , g/m^ EzY is 
divergence free: 

V • a = ^dr{r VGCr) + ^^el^GcJ + + {qt/m, E^) = (5) 

because of variable independence d^^^Ez — dv^^{dz^{r,9, z)) = and dzV^\ = 0. 
Moreover we have vgc — {E x B)/B^, with E = — V<I> and B ~ Bz e^, thus: 

VGC^ = 4- ( ) and vgCb = 4" i^r'^) (6) 
Bz\rJ Bz 

and 

• a = -drir VGC.) + -de{vGC„) = {dr {r {-l/r)de<f) + dg (9,$)) = 0. 
r r r Bz 

(7) 

Therefore, one can write an equivalent conservative equation to the preceding 
Vlasov equation 

dtf + dr{vGC. f) + de{vGCe /) + ^.(^'11 /) + 9.,, (^^Ez = 0. (8) 

This conservative system will be discretized using a conservative semi - La- 
grangian scheme. Following [T] and |16| . we consider two conservative schemes, 
which are fourth order in space: 
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• LAG: LAGrange polynom method, which uses Lagrangian polynoms to 
reconstruct the distribution function. 

• PSM: Parabohc Sphnes Method, which uses cubic sphnes to reconstruct 
the distribution function. 

These schemes are designed to solve conservative models and they allow a di- 
rectional splitting of the Vlasov equation ^ . This equation will be solved by 
using D (dimensions of space) ID conservative steps, discretized by using the 
ID schemes (LAG or PSM). At the continuous level, each ID step has no max- 
imum principle, it is only the solution after all D directional steps, the solution 
of the Vlasov equation, that should satisfy a maximum principle [T] . However, 
high order schemes may create spurious oscillations leading to break this max- 
imum principle. A flux limiting procedure may improve the discrete solution, 
which may be closer to the maximum principle in the sense of showing less 
spurious oscillations. The limiter does not ensure a maximum principle in ID, 
but should decrease the oscillations amplitude created by the scheme by locally 
adding numerical diffusion. That leads us to investigate and compare in details 
the properties of many limiters, which depend mainly on two issues: 

• How to make the high order schemes degenerate into a more diffusive 
scheme? 

• How to detect in the function profile the location where the scheme will 
produce oscillations? 

We first describe a Hermite formalism proposed by [3] applied to the PSM and 
LAG schemes leading to a finite volume form equivalent to the original semi- 
Lagrangian schemes. This formalism is very efficient to introduce limiters in the 
PSM or LAG schemes. Consequently, we will compare some limiters, focusing in 
particular on the OScillations limiter (OSL) proposed by [4]. We also propose 
and investigate a new limiter (Slope Limited Spline, SLS), based on classical 
slope limiting methods. We evaluate the performances of each limiter using a 
benchmark developed in the Gysela code, which runs a 4D drift-kinetic model 

m- 

The outline of this paper is the following : in section [5] will be recalled some 
important properties concerning the conservative form of the Vlasov equations. 
Then the Hermite formalism applied to LAG and PSM schemes will be explored. 
In section 21 some limiters are described and they are further investigated in 
section [S] in the context of a 4D drift-kinetic model. At last we will comment 
on numerical results. 

2 Numerical schemes for the Vlasov equation 

2.1 Directional splitting of the advection problem 

In a phase space of dimension £), we consider a distribution / which is advected 
by a velocity field a. The model taken into account satisfies V ■ a = 0. 

dtf + V, • (a/) = 
V ■ a = (9) 

f{x,t)>0 

INRIA 
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For instance in cylindrical geometry, x = (r, 0, considering the 4D 

problem dealt by the Gysela code. In the next sections, we use a directional 
splitting following yy by solving the conservative system ([3 by D separate ID 
problems for each phase space direction which are still under a conservative form. 
So formally, we will consider the problem (|10p for each of the D directions. The 
generic direction is named x. 

r dia{x,t)f{x,t)) _ 

<GM+,xeR,a(x,t) dx (10) 

In this context, we don't have in general Vfc £ = 0, but only 

V-a = 0. 

2.2 Distribution function and phase space 

We divide one direction of the phase space, generically x with a constant step 
Ax to get a regular mesh. The cells are numbered by an integer i from to N 
and the cell faces by an one-half integer i ± 1/2 (see fig. [1]). Hence we have N+1 
cells and N+2 faces. 

X,_3/2 2;j-l/2 Xi+i/2 X,+3/2 

I ' — rr^ — ; — ^ — m I ^' 

-1/2 1/2 3/2 N-1/2 N+1/2 

' i~o ' i~i ' i=N-i ' ~ ' 

Figure 1: Mesh grid on the x phase space direction (top) and beginning (left 
bottom) and end (right bottom) of the mesh grid. 

We note the distribution function at the time t and at the position x in the 
phase space: g(t, x). We discretize the time space with a constant time step At. 
Writing t" — nAt, we then note g"{x) — g^f^^x) the value of the distribution 
function at t". Using the previous discretization, we define the distribution 
function at the cells faces as (?"_)_]^/2 ^ 5'"(2^i+i/2) (fig- 12. 2p . At last, we define 
the average of the distribution function on one cell i at t" by: 



A, 



1 rXi+i/2 

5"(x)da;, i = {)...N 

3=1-1/2 



.9i-3/2 5i-l/2 5i+l/2 5i+3/2 

9i-i 9i 9i+i X 

Figure 2: Normative example for the distribution g in the cells and at the nodes 
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2.3 Conservative semi-Lagrangian scheme principle 

The mass conservation in a lagrangian volume of the phase space between t" 
and reads as follows: 

gix,r+^)dx^ [ g{x,t")dx (11) 
yoi"+i Jvoi" 

with Vol" = {X(a;"+\r)|X(x"+\r+i) £ Vor+^} 

where X(a;,i"+^) describes the characteristic curve which passes by — 
X(x"+\t"+i) at Thus X(a;"+\r) is the point by which the trajectory 

passes at t" such as this trajectory also passes by a;"+^ at The character- 

istic curves are obtained by solving the following equation: 

dX(x,t) , , 



by the cells faces a;i_|_i/2 at t"+^ which are noted: Xi+xji = X(a;"^j^^2) ^"~''^) ^-nd 
we introduce x^^^j^ which is the point on each characteristic curve at the time 
t": 

The conservation equation (|lip can thus be written using a ID discretized 



form: 



where 



Ax 5?+' = / 9^^\v)dy = [ '^^'^ g"{y)dy (12) 

-1/2 



2^j+l/2 ~ 3:^4-1/2 



= Vor+' = Ax 



In the conservative semi-Lagrangian formalism, the right hand side of equation 
is numerically computed as follows: 



Ax gr+i = / g"(y)d2; - G(x*+i/2) - Gix*^,,^) (13) 



-1/2 



where G(a:;) is the cumulative or primitive function of g defined as: 



G{x) = r 



9{y)dy- 

/2 



This primitive function can be computed exactly at each cell face of the mesh: 

i 

G(x,+i/2)=G(x_i/2) + E^^ ~9k- 

These values at faces are then interpolated by basis functions to obtain an ap- 
proximate reconstruction Gh{x) of G{x) for any a;. 

For instance, the PSM scheme uses cubic splines and the LAG scheme uses 
third order lagrangian polynoma as interpolation functions to obtain the recon- 
structed function Gh{x). 
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2.4 Finite volume form equivalence 

The equation ([T^ can be split in three terms: 

r^i-i/2 rXi+i/2 rx'-i/2 

Ax-g^+' = / 9'\v)dy+ 9"{y)dy+ g^y)dy{U) 

^*-l/2 



<Pi-l/2 Qi^^ -4>i+l/2 

We name 4>i+i/2 tie following quantity: 

'Pr+i/2 = / 9"{y)dy. 

We call it 'flux' since it represents the algebraic quantity which is carried through 
the node by identification with the finite volume formalism. 

The equation is a finite volume equation (see fig. [3]) meaning that the 
new value of the distribution function g^^^ in the cell i is the sum of its value 
at time and the incoming or outgoing flux 4>i±i/2- 

-g^^^ = -g^-{tl±lll_^lzlllY (15) 





1 

1 1 


Xi-l/2 Xi+i/2 






1 1 


^ ,y 

^1-1/2 ^1+1/2 

1 K\\\N\\\\) 


1 1 

1 1 


\ X ' 






1 







0.-1/2,,-— 1 ff? ^^^^^^ 

4>, + \/2.j+=i 

Figure 3: Conservative evolution of a finite volume scheme 

As for the semi-Lagrangian formalism , the integrals to be computed to 
obtain the fluxes can be approximated by a reconstruction of the primitive of 
the distribution function Gt'- 

rxi+1/2 

4+1/2=/ g"iy)dy Ghix^+i/2) ~ Ghix*^i/2)- (16) 
2.5 Time scheme used in the Gysela code 

The time scheme has been modified, because the classical second order leap- 
frog algorithm used in the Gysela code is not robust enough to deal with the 
conservative semi-Lagrangian schemes considered here: 

-n+l ^ -n-l _ f C+1/2 ~ C-1/2 
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Moreover, this time scheme enforces a constant time step because it involves 
three different time steps, which is quite restrictive for high iteration numbers 
simulations. We thus turn the time scheme to a Predictor- Corrector (or Runge 
Kutta order 2) method, which allows to use a variable time step At" computed 
at time t". It is computed according to a CFL like condition necessary for the 
finite volume scheme stability which assesses that the maximum displacement 
in the domain is less than a fraction of the cells size (using a regular mesh): 

Ar = CFL mm { ^^^) (17) 

d=i,D ymax(a2(a;)) J 

with Axd the space step in direction d and a^' the velocity at time in the 
space direction d £ [1, CFL is a coefficient which < CFL < 1. 

Remark 1. The finite volumes scheme form psp and the semi-Lagrangian 
scheme (jl3p are strictly equivalent, since the displacement is restricted to CFL < 
1. The finite volumes scheme (jlSp is not defined for displacements bigger than 
one cell, with CFL > 1, although the semi-Lagrangian scheme p3p could be 
written for any time step. However, the stability of both schemes with CFL > 1 
in a general situation is not demonstrated. 



Predictor- Corrector Algorithm: 

• At beginning of the iteration at t", we compute At" according to ([TT]). 

• Prediction step : we compute a order 1 in time approximation of the 
solution at time t"+i/2 with half a time step according to values ^" and 
cell faces fiuxes (j)^ i+1/2 direction d G [l,D] at time t": 



1+1/2 _ _ ^d,i+l/2 ^dA-1/2 '' 

Axd 



d=l 



We compute the electric potential at same time t"+^/^: 

^ -Vi • (noVi$"+l/2) + -^($"+l/2_ < $n+l/2 = / ,g«+l/2rf^rf„|| 



Buj kT„ 



Correction step : we compute an order 2 in time approximation of the 
solution at time t"^-^ according to values at time and cell faces 
fluxes 't'^d^l_{% in all direction d G [1,-D] at time 



9^' = 5? - E 



rf=i 



rn+1/2 _ J,"+l/2 
^d,i+l/2 'Pd,i-l/2 

Axd 



n+l. 



We compute the electric potential at same time f 
-^Vi • (noVi$"+i) + -^($"+1- < $"+1 >0) = / r^^d^idv\\ - no 

BiOi Kip -r J II 
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3 The Hermite formalism applied to PSM and 
LAG 

Two schemes can be used to compute the new distribution function (reconstruct- 
ing the flux or the primitive function): 

• The Parabohc SpUnes Method called PSM, which uses cubic splines func- 
tions for the interpolations. 

• The Lagrangian method (LAG), which uses third order Lagrange poly- 
noms for the interpolations. 

The PSM scheme [IJ or the LAG scheme |4j are conservative semi-Lagrangian 
schemes that only differs by the interpolation functions used for the reconstruc- 
tion step. We aim to study some limiters for both schemes and using the Hermite 
formalism. The Hermite formalism is a generic formulation for the interpolation 
polynoms. Indeed, the scheme LAG or PSM in this formalism are set only by 
the way of computing the distribution function at the faces 3^+1/2- Thus we can 
easily use them simultaneously in a code. 

We rewrite hereafter this generic formulation, based on the Hermite formal- 
ism for the conservative schemes PSM and LAG, proposed in fi]. First, we give 
the expression of the flux 0i+i/2, then we propose an application to the PSM 
and LAG schemes. 



3.1 Computation of the flux 0i+i/2 with the Hermite for- 
maHsm 

Following |16j and [1] , we reconstruct the distribution function g with a second 
order polynom Pj-, which interpolates the distribution function in the cell k. 
We first assume that the reconstructed function is continuous at the cell faces. 
We name the value of g at the cell faces gk+i/2- At the end of the section, 
we will give the general formula for a discontinuous function at the cell faces. 
The distribution function is approximated by a second order polynom in cell fc, 
which corresponds to interpolate the primitive of the distribution function with 
a third order polynom: 

Vfc, Pk{z) = c^^^ + b^^^h. + a^^h"^ with z G [0, Ax] 

The distribution function g is continuous at the faces. That implies: 

{Pk{Q)=9k-l/2 

|Pfc(Aa;) = gk+i/2 

A third condition on the polynom comes from the mass conservation in the 
cell k, that reads: 



-7— / Pk{y)dy = gk 



Therefore, the polynomial coefficients can be written: 



a 



^''^ = (35fe-i/2 + 35fc+i/2 - 6gk) /Ax^ 



fcC^) - {-Agk-1/2 - 25fe+i/2 + 6.gfe) /Ax (18) 

9k-l/2 



RR 11° 7467 



12 



Guterl et al 



Finally, Vx G [a;fe-i/2, a^fc+1/2]! Pk{x) is an approximation of the distribution 
function g in the cell k. This approximation depends on .9fe-i/2, 9k+i/2i a-i^d g^. 
Using this interpolation, we are able to evaluate the flux (t>k+i/2- 



yfc+1/2 



Xk+1/2 

5"(2/)dy. (19) 



We take the additional assumption that the displacement at the cell face 
k + 1/2, i.e. ak+1/2 = Xk+1/2 - x*k+i/2 satisfies: 



|afc+i/2 



< Aa 



It is necessary to know in which cell j {j = k or j = k + l since the displacement 
l'^/c+i/2l ^ ^2;) is located the foot x'^^-^^^^ of the characteristic that passes by 
Xk+1/2 at time t^^^. We can thus give an approximation of the fluxes (fTO|) 
towards the face j + 1/2 based on the set of polynoms Pk : 

3i I xl^if2 G [2^^-1/2) 2:^ + 1/2] 

(/),+i/2= I g'\Y)dY^ I Pj{Y)dY 

•^^fc + l/2~^j-l/2 ■^^fc + l/2~^3-l/2 

with the change of variable Y — y — 6 [0, Ax]. 

We exhibit hereafter the flux 4>k+i/2 function of the gk-i/2^ 9k+i/2i and gk. 
We note jfe+1/2 the cell where is located the foot of the characteristic passing 
by the face k + 1/2: 

The index jk+1/2 indicates from which cell the flow is coming. For instance, a 
negative displacement < on the face /e+ 1/2 means that the flow comes 

from the cell fc + 1, so that ik+1/2 = fc + 1- 
So we have: 

• If ^k+i/2 < Xk+1/2 then 

- ifc+1/2 = k 

4'k+l/2 — Jxl^-^^.-^~Xk+i/2 

• If ^fc+1/2 > Xk+i/2 then 

- jk+1/2 = k + l 

■^fc+1/2 ■^fc+1/2 

We introduce 

. _ Xu+\I2 - 2^Jfc+i/2-l/2 

Ax 

thus 5 — Q ox 8 — \ indicates the upwinding direction. 
By introducing a normalized displacement 

_ 2^fc+l/2 ~ 2^fc+l/2 
^~ A^ ' 



Ja;* , 2:1,4-1 /2 k\ ) 
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we can write the flux: 



J{S-l3)Ax 



That leads to 

(t>k+l/2 



Ax 

with 



''Jfc+l/2 2 

^93fc+i/2-V2+3gj,,^^^3+i/2-6gj,^_^^^3 



^ ''ifc+1/2 



By ordering differently the polynomial expression ((20)) and replacing the 
coefficients by their values, we try to get an equivalent formulation to the flux 
expressions proposed by [1]. 



1/2 



5,.,v2+i/2 (/35 + /3'(l-3<5)+/33) 
5,.,i/2 (/3'(-3 + 6<5)+/33(-2))] 



with 



/3 = ^^il^^^^^G[-i,i] 



(5 



if a;fc+i/2 < 2^^+1/2 

1 if a;/c+i/2 > ^*k+\i2- 



It's obvious that /3 depends on k hence formally ji — Pk+i/2 

Positive displacement a > i.e. 5 = 1 

We have jk+1/2 = hence 

'/'fc+l/2(/3) . 2/0 

- .9fc-l/2 [P [P - l)j 

+ 9k+l/2 (Pil - P)^) 

+ gk{P'{S-2P)) 

Negative displacement a < i.e. 6 = 

We have jfc+1/2 = fc + 1, hence 

0fc+l/2(^) , 1\2\ 
= .9fe+l/2 [PiP + 1) ) 

+ 5fc+3/2(/3'(l+/?)) 

+ .9fe+3/2(/?'(-3-2/?)) 

Thus, we have established a generic formulation for the flux considering 
displacements smaller than one cell. We have assumed that the distribution 
function is continuous at the cell face. This generic expression is equivalent to 
the formulation proposed by 
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Hermite formalism for a discontinuous reconstruction 

Performing the same calculations, we can extend the previous formalism to a 
distribution function which is not continuous at the faces. We name .gj|!^2/2 ^^^'^ 
9k+i/2 respectively the left and right values of the distribution function at the 
cell face k + 1/2 (see fig. H]). In general that means: 

7^ -9^+1/2 

The conditions on the polynoms at the cell faces are thus changed to: 

4-1/2 Vl/2 

4 — 




Figure 4: Asymmetric node 

\Pfe(Ax)=g+^,/2 
That leads to the final formulation: 

+ .W/2(/5'M + 65) + /33(-2))] (21) 



where 



_ Qi+l/2 ^ Xk+l/2 - ^u+xjl , , 

^~ Ax ^ Ax ^ ^ ^' 

S=l ^ if •^fc+i/2 < a;*+i/2 

1 1 if > 4+1/2 



Positive displacement a > i.e. S — 1 

We have jk+1/2 = k, hence 

0fe+l/2(/3) 



5fe-i/2(/3'(/3-l)) 
+ 5fc(/^'(3-2/3)) 



Negative displacement a < i.e. 6 = 

We have jk+1/2 = + 1, hence 



+ .9,++3/2 + /?)) 

+ 5fc+i(/32(-3-2/3)) 
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3.2 PSM scheme with the Hermite formahsm 

3.2.1 Hermite formalism using splines interpolation 

The type of scheme expressed in the hermite formahsm (PSM or LAG) only 
depends on the manner the face values g^j^i^2 ^^'^ 9k+i/2 computed. Here 
the PSM scheme has the property that the reconstructed function has continuous 
derivatives of the distribution function g{x,t) at cell faces [TB]. Remembering 
that the distribution function g is reconstructed by a second order polynom 
g(x) « Pi{x) in the cell i, the continuity of the derivative at face i + 1/2 can be 
written: 

li-^A.- |y=o 
Using the polynom coefficients expression ([T5| we obtain: 

dP,{Y) . dP,+i{Y) , 

<^ 2aiAx + bi = 
^ 3i-i/2 +4gi+i/2 +3^+3/2 H9i+9i+i) 

This PSM formulation regardless of the boundary conditions is equivalent to 
the semi-Lagrangian PSM formalism, used in [1 1 for instance. A rigorous proof 
of the equivalence between the two formulations is furnished by [4J, except for 
the boundary conditions. 

3.2.2 Periodic boundary conditions 

In this section, we present the way of imposing the boundary conditions. 
Extending boundary conditions for the PSM scheme to the Hermite formalism 
is not simple, especially to get a complete equivalence between the Hermite 
formalism and the semi-Lagrangian formalism (|13p [T] . 
We named G the primitive function which is defined as: 

G{x) = f g{x)dx (22) 

We also define the mesh fitted to a periodic domain(fig. [5]): 

-1/2 = 1/21/2 3/2 N-l/ff + 1/2 = -1/2 
I 1 1 -. 1 1 1 . ^ ^ 



1 N N+l= 

J 



2tt 

Figure 5: Mesh for a periodic domain 



The periodic boundary conditions for the semi-Lagrangian scheme are ap- 
proximated by conditions on the primitive function derivatives: 

G'(a;_i/2) = G'{xN+i/2) 

G"(x_i/2) = G"ixN+l/2) 
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which is equivalent to set continuity of the distribution function and its first 
derivative: 

\fi''(a;-i/2) = g'{xN+x/2) 

The periodic boundary conditions for the Hermite formahsm are then obtained 
by setting the same constraint on the polynomial reconstruction: 

/Po(0)=Pjv(Aa;) 

\ dPo(Y) I _ dPjv(r) I 

^ \ 9-1/2 = 9N+1/2 

\-ig-i/2 - 2c/i/2 + 6go = 2gN-i/2 + ^9n+i/2 - ^Sn 

The complete linear system to solve to obtain values at nodes is then the fol- 
lowing 

9i-i/2 + ^9i+i/2 + 5i+3/2 = 3(5i + i = Q---N -1 

9-1/2 = 5^+1/2 

^-4fi(_i/2 - 2.91/2 + 650 = 2ffAr_i/2 + 45a,+i/2 - Qqn 

gi-i/2 + 45i+i/2 +gi+3/2 = 3(5i + gt+i),i = Q---N -2 
=^ i 9N-3/2 + ^9n-i/2 + 9-1/2 = 3(5Af-i + 9n) 
^9n-i/2 + 4g_i/2 + 9i/2 = 3(50 + 9i) 

The matricial system [A\X = B of dimension TV + 1 is finally: 



/4 1 
1 4 



Vl 





1 







l\ ( 9-1/2 \ 
9l/2 



14 1. 

Q ' 'I. W-1/2J 



( K9N + ,90) \ 

3(30+51) 



3(5Ar-2 +5W-l) 

V 3(5Ar_i + 5jv) j 



3.2.3 Non-Periodic boundary conditions: natural conditions 

The natural boundary condition for the PSM scheme is approximated by vanish- 
ing the second derivative of the primitive function, what corresponds to annulate 
the first derivative of the distribution function 5, at the boundaries: 



G"(x_i/2) = 
G"(xjv+i/2) = 



The corresponding natural boundary conditions for the Hermite formalism are: 

-45-1/2 - 251/2 + 650 = 



dPo(Y) I _ „ 

dPNjY) I _ f^ 



25JV-1/2 + 45JV+1/2 - 65JV = 
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The system to solve to obtain the values of the function at the nodes is then: 



gi-i/2 +4gi+i/2 +gi+3/2 = 3(.gi + = 0---A^- 1 

45-1/2 + 2511/2 = &go 

Considering g-1/2 — go and gN+1/2 — Sn, we solve from the following linear 
system of dimension + 2: 



/4 2 
1 4 



Vo 





1 



0\ / 5-1/2 \ 

.91/2 



gN-1/2 

^' , \ffJV+l/2/ 



( 6(50) \ 
3(go + .gi) 



3(ffjv-i +5^) 
V 6(gjv) y 

V 



3.3 LAG scheme with the Hermite formahsm 

3.3.1 Hermite formalism using Lagrange interpolation 

We outline hereafter the Hermite formulation for the LAG scheme. In order to 
set the distribution function values at the faces g^^-^i2 and g^j^Yi2^ start from 
the Lagrange interpolation of the flux and then, we make some calculations to 
show the equivalence with the Hermite formulation. 

We build the primitive or cumulative function G at the characteristic feet: 



5 



n+l 



G{'A+\I2) ~ G(^j*-l/2) 



G'(a^*+i/2) — 



g'\x)dx. 



-1/2 



The value of G{x*^j^^2) obtained by an interpolation with third-order La- 
grange polynoms of values at cell faces G(xj_|_i/2): 



i+2 



Gix) 



Gj-l/2Lj{x),X ^ [Xi_lf2j ^i+1/2] 

j=i-l 



where Lj are the Lagrange polynoms defined as: 

i+2 



Lj{x) 



n 



X - Xfc„i/2 
2^^-1/2 - a;*;-l/2 



and Gj_i/2 is the value of G at the cell face j — 1/2. 

Some calculations (given in annex sec. 17. ip lead to the expressions ([23| and 
for Gix*^^^^): 
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G«+i/2) = /33(l/6G,_3/2-l/2G,_i/2 + l/2G,+i/2-l/6G,+3/2) 

+ /32(l/2G,_i/2-G,+i/2 + l/2G,+3/2) 

+ /3(-l/6G,_3/2 + G,_i/2 - l/2G,+i/2 - l/3G,+3/2) 

+ Gi+i/2 (23) 

• if ft+i/2 < 0, 

G(a;*+i/2) = /33(l/6G,_i/2-l/2G,+i/2 + l/2G,+3/2-l/6G,+5/2) 
+ (l/2G,_i/2 - G,+i/2 + l/2G,+3/2) 

+ /3(l/3G,_i/2 + l/2G,+i/2 - G,+3/2 + l/6G,+5/2) (24) 

with Gfc+i/2 = G(xfc+i/2) and the function /3i+i/2: 

We omit the subscript for /3 when the context makes the index obvious. 

The values of the primitive function at nodes Gfc_|.i/2 are related to the cell 

centred values of the distribution function by: 

fffc Aa; = Gfc+i/2 - Gk-i/2- 

Thus, by replacing the values Gfc+1/2 function of in such a way only the 
primitive function value Gi_|_i/2 at cell face i + 1/2 remains at the right hand 
side of relations and (j24p , we obtain: 

• if ft+1/2 > 0, 

G(x*+i/2)/Ax = /33(-l/6gr-i + l/3.gr-l/6gr+i) 

+ /32(-l/2gr + l/2gr+i) 

+ /3(l/6gr-i - 5/65," - l/35.'Vi) 

+ Gi+i/2/Ax 

• if A:+i/2 < 

G(x*+i/2)/Ax = /?3(-l/65r + 2/6.gr+i-l/6gr+2) 
+ p\-l/2g- + l/2g-^,) 
+ /3(-l/3.gr - 5/6.gr+i + 1/6.97+2) 
+ G,+i/2/Aa; 

Let us recall the relation (fTC)) between the flux at cell face i + 1/2 and the 
primitive function: 

Therefore, we have an expression of the flux </'i+i/2(/5) function of cell centred 
values {gk)k- 
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if > 0, 

0i+l/2(/3) 

Aa; 



• if < 



/33(l/65r_i-l/3<7r + l/65r+i) 
/32 (1/2.97 - l/25r+i) 
/3(-l/6gr-i + 5/65," + l/3.9r+i) 



'^±^^ = /33(l/65r-2/65r+i + l/6,gr+2) 



+ /32(l/25,r - l/2.gIVi) 

+ /3(l/35r + 5/6gr+i - 1/6.9.'V2) 

Let us remember the generic Hermite formulation for the flux (PT|) function of 
3fc+i/2- 



values at cell faces 



if A+1/2 > then 

'^»+l/2(/3) 

Ax 



/?'(-.9+i/2-2,9r+i/2 + 3<?r) 

/^(5r+i/2) 



if A+1/2 < then 

^i+l/2(^) 

Ax 



/3 (5r+3/2+. 9^1/2 -2.9r+i) 

+ /3'(2.9++i/2+.9r+3/2-35lVi) 
+ /?(5ii/2) 

By identifying in the two last expressions of the flux 0^+1/2 the coeficients 
of each degree of the polynom in variable /?, we obtain a necessary expression 
for values g^_i/2 faces: 

/5fe-i/2 = l/3.9ri + 5/6.9L' - 1/6.9^+1 ^^^^ 



[9k+i/2 - -1/6.91-1 + 5/6ffI? + 1/35?+ 



As a conclusion, we have explicit values of g^_j^^2 ^^'^ 9k+i/2 function of 
cell centred values of the distribution function at time t" with (pS)) , so we can 
directly compute the generic Hermite formulation for the flux ((2T|) for the LAG 
scheme. 



3.3.2 Boundary conditions 

The boundary conditions for the LAG scheme are more simple than the PSM 
ones. Indeed, the coefficient are explicitly depending on the value of the distri- 
bution function g^'. So we can explicitly set values at boundaries. 
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Non-periodic boundary conditions: natural conditions 

'9n+3/2 = V3ffKr+i + 5/Q9N+1 - V65Xr+i 

^ 9N+3/2 = -^/^9n + ^/^9n+i + ^/^9n+i 
' 5ii/2 = l« + 5«-l/65r 

,511/2 = -1« + 5/6flff + 1« 
Periodic boundary conditions 

'.9+^/2 = l/3g]^_i + 5/650" -1/65? 

^ 5:1/2 = -1/653^-2 + 5/6.9Xr-i + 1/350" 
' 5^^.1/2 = 1/355^-1 + 5/65^-1/65? 

.5^-1/2 = -l/65]^_2 + 5/65]^_i + 1/35? 
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4 Limiters 

High-order numerical schenies may create spurious oscillations when stiff profiles 
occur in the distribution function. As usual in the numerical framework of 
finite volume schemes, we employ a flux limiter to lower these oscillations. The 
principle of the limiter is to introduce numerical diffusion when stiff profiles are 
detected, by modifying the flux at the cell faces. We will test limiters found in 
the literature, the so called Oscillation Limiter (OSL) described in and a 
the one proposed in this report that we call Slope Limited Splines (SLS). These 
limiters are mainly provided with the PSM scheme, but some might be used for 
any finite volume scheme. 

4.1 ENTropic flux limiter (ENT) 

The principle of this limiter proposed in [IJ is to make degenerate the fourth 
order PSM scheme to a second order centred flux to reduce the anti-diffusive 
behaviour of the scheme when it occurs. The position where anti-diffusion occurs 
is detected looking at the second order equivalent equation solved by the scheme, 
obtained by a Taylor expansion. This equivalent equation shows a diffusion term 
at second order, which the sign should be positive, leading to numerical diffusion 
and thus stability. When this sign is negative, the scheme is anti-diffusive and it 
is replaced by a centred scheme. The corresponding algorithm is the following: 

Algorithm 

• Computation of the PSM fiux ^^^^/a ^^^^ ^^'^ centred flux ^^^/'j = "^+1/2^^-^^ 

• if (Ci72 - €+1%) i9?+i - 9?) > then the flux (/-f/^^ is supposed to 



be diffusive 

• if (<^f/i72 - Cvs) (^'"+1 - 9?) < then the flux (/.f/*^^ is supposed to 
be anti-diffusive, thus it is switched with the centred flux: 

J,PSM _ ^CEN 
Vi+1/2 — Vi+l/2 

4.2 UMEDA's limiter (UMEDA) 

The 4D advection equation is split in four ID advection equations. Although 
the maximum principle is satisfled by the 4D equation, this principle is not 
fulfilled for each ID equation. Therefore the extrema of the distribution function 
(minimum and maximum) are unknown in ID. The PFC limiter [7J employ 
these extrema, which are not known in the context of a directional splitting. 
Therefore, following Umeda [15], we modify the extrema definition to get a 
limiter working in 4D. The Umeda's limiter have been written with Lagrange 
polynoms. 

Algorithm 
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• We evaluate gminl, gmin2, gmaxl, gmax2 

' gmtnl = max [max {g'^_i,gf) , min {2gf_^^ - gf_2,'2gf - 5^+1)] 
5m*«2 = max [max (gf+i, gf) , min {2gl'_^-^ - g'l^^, 2^^ - g'i_^\ 
gmaxi = min [min , max (2.gf_i - 5^,2, 2gf - gj^^)] 

,5ma2;2 = min [min (gf+i, ) , max (25f_^i - .gf_^2> - fff^i)] 



• We set Qrnin-i Qmax 



I ymin — 

max [0, min (g 

mini 7 gmin2 } 
I gmax — max [gmaxl: gmax2] 



• We define Lf and L, 



l: 



if - 9i > 0, min(2(gj - gmin),gi+i - gi) 
if g-j+i - (ji < 0, max(2(g, - g„,ax), (ji+i - gi) 

I if - > 0,min(2(g,„aa; - - 

I if - < 0,max(2(5„„ - gi),gi - gi-i) 



• We finally redefine the flux: 

G,+i/2{x) = 5, + x{l - x){2 - x){L+/Q) + x(l - + .x)(L-/6) 

We obtain the LAG reconstruction without limiter by setting: 

I -^fe 9k+i - gk 

— 9k — gk-1 

4.3 oscillation Limiter (OSL) 

The oscillation Limiter (OSL) proposed in U is really using the Hermite for- 
malism. It compares the values at call faces g^_i/2 obtained with the LAG and 
the PSM schemes to the value computed with a linear reconstruction. This later 
consists in computing an average at the face of the left and right cells centred 
values. If the values computed with PSM and LAG are not simultaneously up- 
per or lower than the average value, then we take the average value. If not, a 
mixed scheme between PSM, LAG and the average reconstruction is performed. 
The OSL limiter includes a parameter C>1 determining the proportion in the 
average of PSM and LAG fluxes in the limited flux. 

Algorithm 

• Computation of g^_i/2 values for both PSM et LAG schemes. 

• Average value at the k — 1/2 and k + 1/2 nodes: 5™'\y2 ~ (Sfe + 9k-i)/'^ 
and 9r+i/2 = i~9k+9i:+i)/2- 
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• The following formula perform the choice for the limited flux according to 
the regularity of the cell faces values: 

- if i9t-i/2,LAG - 9T-\/2)i9k-i/2,PSM ~ 9^-1/2) > « then 

9k-l/2 ^ 9k~\/2 ^i^^ 

„+ _ „ave 

yfc-1/2 ~ ■y/c-1/2 

+sign{gk-i/2,PSM ~ 9T-1/2) (26) 
min(C|5+_^/2,LAG ^ fffc-l/sl' \9k-i/2,PSM ~ 5^-1/2!) 

- if i9k+i/2.LAG ~ 9t+i/2)i9k+i/2,PSM " 9T+1/2) > t^en 
9k+i/2 ~ 9T+1/2 ^i^^ 

9k+l/2 ~ 9k+l/2 

+sign{gk+i/2,PSM - 9I+1/2) (27) 
min(C|5~_^^/2_^^Q ~ 9t+i/2\^ \9k+i/2,PSM - 9k+\/2\) 

4.4 Slope Limited Splines (SLS) 

The limiting procedure basically aims to cut the oscillations generated by strong 
gradients in the distribution function profile, where high order schemes will 
produce overshoots and spurious oscillations. We propose here to measure these 
gradients and to add diffusion where strong gradients are detected. The diffusion 
is added by mixing the high order scheme with a first order upwind flux. The 
more the gradient is steep, the more we raise the proportion of upwind flux in 
the average with the high order scheme. The evaluation of the gradient is given 
by a function 9 and we estimate the diffusion needed with a function "f{9) G [0, f ] 
based on the minmod like limiter function (fig. [T]): 

Cv2 - 7(e.+i/2) + (1 - 7(e.+i/2)) €rijf 

where 

{lUpwind 9i +9i + i ■ / \9i+i~9i 
<Pi+i/2 = "1+1/2 2^ ~ sign(ai+i/2)—^ 
"1+1/2 = A<a,j+i/2 

We define 9i^i/2 as the classical slope ratio of the distribution which depends 
on the direction of the displacement (fig. 



-9i- 



'i+1/2 



fi+i~^i 

-,r- if /9 



;/,+2-g.+i 



if > 

if a-i+1/2 < 



However, the classical limiters as minmod 7i+i/2 = max{Q, ^'^(^1+1/21 l))i set 7 
to when 9 < Q. That means that the scheme turns to order f when an extrema 
exists, i.e. the slope ratio < 0. These extrema are thus quickly diffused and 
that leads to loose the benefits of a high order method. For SLS, the choice is 
to let the high-order scheme deal with the extrema and only add diffusion when 
high gradients occurs, i.e. the slope ratio 0^0. We also introduce a constant 
K in relation to control the maximum slope allowed without adding diffusion. 
The SLS limiter function is thus to set 7 = 1 for any values of 9, except close 
to \9\ — where strong gradients occurs, see figure [T] 

lz+i/2 = max(O,min(ii:|0,+i/2|,l)) (28) 
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5 Numerical results 

We consider both ID and 4D test cases. The ID advection test case is performed 
with a constant velocity to give a quahtative view on the effects of each hmiter. 
We are particularly interested in advection of profiles with strong gradients, as 
discontinuous initial distribution functions as a step. However, we have to keep 
in mind that this situation of a discontinuous shape advected with a constant 
velocity field does not occur in Vlasov 4D drift-kinetic simulations for many 
reasons: the advections are not constant, the long time constant advection con- 
figuration does not appear, and discontinuous functions (steps) do not exist in 
the Vlasov model. However, it is relevant to investigate advection of discontinu- 
ous functions since Vlasov models leads to stretched structures in the flow that 
makes appear strong gradients, which could be assimilated to discontinuities. 

The 4D drift-kinetic test case is performed with the 5D Gysela code follow- 
ing the benchmark presented in [9]. This benchmark is to simulate instabilities 
growing in the plasma leading to shapes like filaments and vortex. We com- 
pare the results to evaluate the limiters effects on the development of turbulent 
structures. We will focus on two directions, (r, 6*), among the four (r, 0,0, 
since strong gradients of the distribution function appear mainly in (r, 9) planes 
(further details in [1]). We give algorithms used in Gysela to solve the advection 
equation in sec. 17.41 



5.1 Tools 

We propose hereafter some quantitative tools to compare the efficiency of each 
limiter. The norm and the entropy may be used to investigate the diffusivity 
of each limiter, since each should theoretically be kept constant. The total 
energy conservation is related to the treatment of small structures provided by 
the instabilities. In addition, we also introduce the total variation (TV) norm 
to estimate the numerical oscillations created by the scheme. 



Total Variation 

The total variational norm (TV) is used to estimate the rate of oscillations 
produced by a scheme, compare to another. 

• In ID: 

4=0 



In 4D, we restrict the diagnostic to an {r,d) plane: 
TV{g{t)) - / / g{t,r,0)\L2drd9 



1=0 j=0 



''ff"(* + l,j)-.9"(*,j)\' , /.9"(«,J + l)-ff"(*,j)\' 



Ar / V AO 
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norm and entropy 

The norm is defined as: 

L^{f{t)) J {f{t,X,v\\)fdXdv\\, dX = rdrded(j). 

The entropy measures the information created or destroyed by a phenomena. 
We can thus expect that a raising entropy estimates the information lost by 
numerical diffusion. 



S{f{t))^- J fit,X,vii)logifit,X,v\i))dXdv\i, 



dX = rdrded(j) 



Total Energy 

The conservation of the total energy etot writes as the sum of the kinetic 
energy ekm and the potential energy epot defined as follow: 

Etot = Cfcin + epot = I ^m^vlU - feq)dVdv\\ + ^ / ^'t^i'^i ~ no)dX. 



For further details, see [5]. 



Quality factor 

We propose here to define a quality factor to establish a quantitative evalu- 
ation of the limiters efficiency Since the norm should be conserved by the 
Vlasov equation, it furnishes a good estimation of the diffusive behaviour of 
a scheme, because it decreases this norm. On the other hand, the TV norm 
estimates the numerical oscillations, or also the diffusion in a way, but more 
'locally'. The quality of a scheme may be defined as its ability to limit the oscil- 
lations (small TV) with the less numerical diffusion (high norm). Therefore 
we introduce the quality criterion Q defined as: 

= mm 

5.2 Test of the limiters on constant ID advections 

We test the limiters in ID on a step function 0. The step exhibits the problems 
for which the limiters are required. This ID benchmark gives a first overview 
on the limiters capabilities. 

The benchmark consists in solving a constant advection on a periodic domain 
divided in 80 cells. The displacement is set to 0.2 cell per iteration. The 
indicated times are the number of iterations. 



^the distribution function f is set to 1 on a part of the domain else to 
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PSM,t=10 PSM, t=200 PSM,t=1000 




10 20 30 40 so 60 10 SO 30 40 50 BO 10 20 30 40 50 60 



XXX 

Figure 8: Constant advection on a step with the PSM scheme and the PSM 
scheme with entropic hmiter ENT. The domain is meshed on 80 ceUs with 
periodic boundary conditions and the displacement is set to 0.2 ceU per iteration. 

PSM-ENT We observe on figure[5]that the ENT hmiter cuts off the spurious 
oscihations of the PSM scheme at the left side of the discontinuity. However, at 
the right side, oscillations are weakly damped. 
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Figure 9: Constant advection on a step with the LAG scheme and the LAG 
scheme with the Umeda hmiter. The domain is meshed on 80 cehs with periodic 
boundary conditions and the displacement is set to 0.2 cell per iteration. 



LAG-UMEDA We observe on figure [5] that the LAG scheme is much more 
diffusive than the PSM scheme without really better satisfying a maximum 
principle. The under/overshoots of the LAG scheme are cut off by the Umeda 
limiter. Although the LAG scheme with Umeda limiter respects a maximum 
principle, the step signal is quickly diffused and it seems not appropriate for 
long time simulations. 
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Figure 10: Constant advection on a step with the PSM scheme and the PSM 
scheme with the OSL limiter with different values of parameter C. The domain 
is meshed on 80 cells with periodic boundary conditions and the displacement 
is set to 0.2 cell per iteration. 



PSM-OSL We observe on figure [TU] that the OSL limiter reduces the oscilla- 
tions without introducing much diffusion. However, there is still under/overshoots 
at both sides of the discontinuity and a slight offset of the solution occurs com- 
pare to the exact solution, but it seems not increase with time. The constant C 
of the OSL limiter has almost no influence on the results. 
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Figure 11: Constant advection on a step with the PSM scheme and the PSM 
scheme with the SLS hmiter with different values of parameter K. The domain 
is meshed on 80 cells with periodic boundary conditions and the displacement 
is set to 0.2 cell per iteration. 



PSM-SLS We observe on figure [TT] that the SLS limiter cuts off oscillations 
at the right side of the discontinuity and keeps the overshoot at left side. The 
results of SLS limiter with K — 1 is too much diffusive (equivalent to a minmod 
limiter). With K — 5, the result is almost the same as with = 10 and the 
accuracy at discontinuities (slope of the solution at discontinuities) is almost 
the same as with the PSM scheme. 
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Figure 12: Constant advection on a step with the PSM scheme, the PSM scheme 
with the ENT hmiter, the SLS hmiter with parameter K — 5 and the OSL 
hmiter with parameter C = 2. The domain is meshed on 80 cells with periodic 
boundary conditions and the displacement is set to 0.2 cell per iteration. 



Comparison of all limiters for PSM It is interesting to observe the action 
of each different limiter for PSM on the oscillations. The ENT limiter cuts off 
the oscillations at left side of the discontinuities, the SLS K — 5 limiter cuts off 
the oscillations at right side of the discontinuities and the OSL limiter reduces 
the oscillations at both sides. 
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t t 

Figure 13: norm (left) and total variation (right) as a function of the time 
for the different methods applied to a step with ID constant advection on a 
periodic domain. 



16 




OSL, C.5 



Figure 14: Quality factor for all the methods as a function of the time for the 
different methods applied to a step with ID constant advection on a periodic 
domain. 

Quality factor The quality factor(fig. [H)) . as defined previously in as 
the ratio of norm over TV norm, shows three groups of curves: 

• On the lower part, the PSM scheme curve shows that even conserving well 
the norm, the TV norm is much higher than the other schemes because 
of spurious oscillations. The quality factor is then poor for PSM. 

• The group of four curves at the middle (LAG, PSM-ENT, PSM OSM 
with C=2 and C=5) are conserving the norm with the same order (best 
PSM-OSL C=2, worst LAG) and with almost the same level of oscillations 
according to the TV norm. The quality factor is thus equivalent for this 
group of schemes. 

• On the upper part, LAG-UMEDA and PSM-SLS K=l (minmod limiter) 
are very diffusive according to the norm, so they kill all the oscillations 
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and the TV norm is very weak. The quahty factor is high because these 
diffusive schemes have too low TV norms and cannot properly be compared 
with the others. However, the PSM-SLS K=5 scheme conserves the 
norm in a comparable way with the other schemes, but with a very low 
TV norm. The quality factor is thus good for this scheme and can be 
compared with to the other schemes. 

As a conclusion for this comparison of all schemes on a linear advection of 
a step function test case, we could say that the LAG-UMEDA and PSM-SLS 
K=l (minmod limiter) schemes are too diffusive to be compared with the other 
schemes considering the quantitative quality factor (|29p . However, this way of 
quantitative comparison seems to fit well for this test case with the qualitative 
comparison or the "feeling" we might have by observing the plotted results in 
Figure [TH 
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5.3 4D Drift-kinetic model 

In this section, we evaluate the hmiter capacities on a 4D drift-kinetic model 
simulation of instabilities with the Gysela code [9J. The benchmark providing 
instabilities is described in [T]. We simulate this test case with two mesh res- 
olutions: a low resolution to provide results about all the limiters and a high 
resolution to get refined results with the best schemes to test also robustness. 
However, let us first do a quick review of the 4D drift-kinetic model and the 
time scheme used for high resolution simulations. 

We recall hereafter the model described in [8j . The geometrical assumptions 
of this model for ion plasma turbulence are a cylindrical geometry with coordi- 
nates (r, 6, z, U||) and a constant magnetic field B = e^, where is the unit 
vector in z direction. In this collisionless plasma, the trajectories are governed 
by the guiding center (GC) trajectories: 

dr de dz q, 

with vgc = {E X B)/B^ and E -V$ with $ the electric potential. 

The Vlasov equation governing this system, where the ion distribution function 

is /(r, 0, z, i), is the following: 

dtf + VGcAf + VGCedef + v\\dj + ^E,d,^^ / = 0. (31) 

TTi'l 

This equation is coupled with a quasi-neutrality equation for the electric poten- 
tial <i>(r, 6*, z) that reads: 

with rii = I f{r, 6, z, W||)df y and constant in time physical parameters no, ^o, 
Te and e. 

Let us notice that the 4D velocity field a — (wgc^ , wcCe 1 j "z/"!; E^Y is diver- 
gence free: 

V • a ^drir VGCr) + -^deivcce) + dzV\\ + a^,, {q/nii E^) = 
because of variable independence dy^^E^ = d^^^ {dz^{r^ 9, z)) — 0, dzV\^ — and 

VGCr = TT^S^ ^^^"^ VGCa/r = —^dr<^, 

r Bz r Bz 

such that ^ ^ 

-drir VGCr) + -de{vGCe) = 0. 
r r 

Therefore, one can write an equivalent conservative equation to the preceding 
Vlasov equation ([^ : 

0. 



dtf + drivGCr f) + deivGCs /) + 5.(«|| /) + a,,,, (^^Ez f 
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5.4 Low resolution simulations 

The low-resolution simulations run on a 64x128x16x16 grid (iV^ xNgX x Ny^^ ) . 
The resolution is too low to consider the total energy as a relevant measure. 
The entropy, norm and TV norm are used to gauge the limiter's effects. We 
present in Figures (TSllini HZl three states of the instabilities: the linear phase 
(left), the beginning of the non-linear phase (center) and the strong non-linear 
phase (right). 

In figure [151 '^^ see the difference of behaviour of the PSM scheme based 
on a spline reconstruction method and the LAG scheme based on Lagrangian 
polynoms. The LAG scheme diffuses much more the solution even at early 
times, and this is accentuated using the LAG-UMEDA scheme which makes 
disappear even the big structures. However, the PSM scheme solution shows a 
lot of small structures, actually at the scale of the mesh, which are suspected 
of coming from spurious oscillations already observed on the ID step test case. 
Moreover, these oscillations often violate the theoretical maximum principle of 
the solution and may lead to crash the simulation. 

In figure [TCI we see the effect of the ENT limiter on the PSM scheme. It 
reduces a lot the oscillations with a figure showing much less small structures, 
but showing a more diffusive behaviour. The PSM-SLS K=l scheme is extremely 
diffusive and makes disappear even the big structures. 

In figure [T71 we compare the PSM scheme result with PSM-OSL and PSM-SLS 
limiters results. Both limiters reduces efficiently the oscillations. The limiter 
PSM-SLS seems a little more diffusive than PSM-OSL. 
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Figure 16: Solution in a (r, 6) plane with NrXNg = 64 x 128 cells (low resolution) 

of a 4D test case with PSM. PSM-ENTand PSM K=l methods. 
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Figure 17: Solution in a (r, 6) plane with NrXNg = 64 x 128 cells (low resolution) 

of a 4D test case with PSM. PSM-OSL C=2 and PSM-SLS K=5 methods. 
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Figure 18: Quality factor (left top), entropy (right top), norm (left bottom), 
TV norm (right bottom) for the 4D simulation low resolution 64 x 128 x 16 x 16 
with all the methods presented. 



In fLgure[TS]we see the effect of each limiter on the quality factor, entropy, 
norm and TV norm. 

• The LAG scheme gives very poor results in the conservation of the 
norm, but shows a medium quality factor. 

• The LAG-UMEDA and PSM-SLS K=l schemes are stih diminishing the 
TV norm at a very low rate as for the ID step test case, because they are 
very diffusive. Even if the quality factor is good and the norm close 
to the other schemes, the qualitative results of figure [T51 and \W\ show that 
even big structures are not well captured with these schemes. 

• The PSM scheme without limiter has still a very different behaviour com- 
pare with the others : the norm and entropy plots show that this 
scheme is less diffusive and the TV norm plots confirms the qualitative 
2D profiles [T7] that it is much more oscillating than the other schemes. 
Thus the quality factor is lower than the other schemes. 
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• The PSM scheme with hmiters ENT, OSL, and SLS K=5 show equivalent 
medium quality factors and TV norms and conservation of entropy and 
norm for this coarse benchmark. 

Conclusion of the 4D results with low resolution section 

The results obtained for this 4D test cas low resolution confirms those ob- 
tained for the ID linear step test case. The PSM scheme shows oscillations as 
the TV norm is high compare to the other schemes, but preserves the better 
norms and the entropy. The PSM-SLS K=l and LAG-UMEDA schemes 
are diffusive and seem experimentally of lower order of accuracy than the other 
schemes. The PSM-SLS K=l limiter, PSM-OSL limiter and ENT hmiter are 
showing a similar level of accuracy. The high resolution simulation results fol- 
lowing will permit to go further in these investigations. 
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5.5 High resolution simulation 

We propose in this section to compare PSM, OSL and SLS with an high reso- 
lution mesh relatively to the previous resolution. The mesh is constituted with 
NrXNgX X iVi,|| = 256 x 512 x 32 x 16. We consider that this mesh is refined 
enough such that conservation of the total energy becomes relevant to estimate 
the limiter quality, what was not the case for the low resolution mesh used in 
the preceding section [531 

We have chosen standard values for PSM limiters according to the preceding 
results : C=2 for the OSL limiter which depends weakly on this value and K=5 
fot the SLS limiter which seems to be the minimum value to limit diffusion (see 
results with K=l in section [^3]). 

The results are presented in three ways : full (r, 9) cut planes, zoom in smaller 
boxes in these cut planes to better see the infiuence of the limiters on small 
structures and ID diagnostics with the total energy, the quality factor (pS)) . the 
entropy, the and TV norms. 
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Comparison of results in (r, 9) cut planes 

In figure [inland [201 we present at different times (r, 9) profiles of the distribu- 
tion function in the non-linear phase. The influence of the SLS or OSL limiters 
is weak at the beginning of the linear phase, i.e. time t — 2290. Afterward, a lot 
of small structures as filaments develop, where steep gradients exist. Therefore 
for later times, the differences between schemes results increase because of the 
different actions of SLS and OSL limiters. 

• The PSM scheme shows spurious oscillations very quickly in the non- 
linear phase. As seen in the ID step test-case and in the 4D test-case with 
low resolution, the PSM scheme develops oscillations when transporting 
discontinuities or steep gradients which occurs in this fiow where filaments 
and vortex develop. Unfortunately, these oscillations lead to the crash of 
the PSM simulation at final time presented t = 3637. 

• The PSM-SLS scheme does almost develop no oscillations even in the 
latest time. The price to pay is a qualitatively more diffused profile of 
distribution function. 

• The PSM-OSL scheme does neither develop oscillations at early times of 
the non-linear phase. The OSL limiter produces less diffusion than the 
SLS limiter. However, at late times the same kind of oscillations than the 
PSM scheme develops. 
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PSM, t=2290 



SLS, K=5, t=2290 



OSL, C=2, t=2290 






PSM, t=2885 



SLS, K=5, t=2885 





OSL, C=2, t=2885 




Figure 19: 4D simulations on Nr x Ng — 256 x 512 cells (high resolution) 
along (r, 0) with PSM (left), (center) SLS, K=5 and OSL (C=2) (right) at three 
different times: t=2290 (top), t=2885 (bottom). 
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PSM, t=3115 SLS, K=5, t=3115 OSL, C=2, t=311 5 




PSM, t=3637 SLS, K=5, t=3637 OSL, C=2, t=3637 
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Figure 20: 4D simulations on Nr x Ng = 256 x 512 cells (high resolution) along 
{r,9) with PSM (left), SLS, K=5(center) and OSL, C=2(right) at three different 
times: t= 3115(top), t= 3637 (bottom). 
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PSM, t=3637 SLS, K=5, t=3637 OSL, C=2, t=3637 




Figure 21: Zoom on the zone A of the high-resolution (r, 9) profiles (fig. [TO]) 
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Figure 22: Zoom on a the zone B of the high-resolution (r, 6) profiles (fig. [T^ 
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Comparison of ID integrated quantities We present in figure[53]the time 
evolution of the total energy, the entropy, the and TV norms and the quality 
factor (Uni)- The results are for the PSM scheme with no hmiter, the PSM-SLS 
scheme with K=5 and K=10 and the PSM-OSL scheme with C=l and C=2. 

• The PSM scheme conserves the total energy the better until the simulation 
crashes down. Entropy and norm show that this scheme is the less 
diffusive but with the highest TV norm, which comes from the spurious 
oscillations seen on (r, 9) profiles and leads to a bad quality factor. 

• On the contrary of what expected looking at (r, 9) profiles, the entropy 
and norm show that OSL limiter leads to a little more dissipation than 
the SLS limiter. The total energy is a little better conserved with the 
SLS limiter than with the OSL limiter at the beginning of the simulation, 
but it constantly decreases while the total energy obtained with the OSL 
limiter stabilise at late times. 

• The quality factor is a httle better with PSM-SLS than with PSM-OSL. 
However, both achieve the objective of stabilising the PSM scheme oscil- 
lations. 
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Figure 23: Total energy (left top), entropy (right top), norm L? (left bottom) 
and TV norm (right bottom) for the 4D simulations (high resolution), quality 
factor Q on the lower part. 
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6 Conclusion 

In this paper, we have investigated conservative schemes for a Vlasov 4D drift- 
kinetic model for plasma turbulence. These schemes should be accurate to bear 
long time simulations, so we investigate fourth order schemes. However, high 
order schemes usually experience difhculties when dealing with steep gradients 
in the flow, which happens when turbulence develops thin structures or vortices. 
In particular, spurious oscillations may appear which are not damped by Vlasov 
models. 

We thus have tested limiters to cut off these oscillations. It consists in adding 
diffusion to the high order scheme at locations where oscillations may appear, 
i.e. at steep gradients. Then two questions show up: how to introduce diffusion 
in the scheme and how to detect the location where diffusion is needed. Of 
course, it should be done efhciently without loosing the accuracy of the scheme 
by introducing too much diffusion. We have tested the LAG and PSM schemes 
and the Umeda, ENT, OSL, SLS limiters, which use quite different strategies. 
We have tested all the schemes on a ID step linear advection test-case and on 
a 4D drift-kinetic model test-case. 

Another question is how to quantify the action of limiters and to propose a 
quantitative way of comparison among them. In a classical way, we have checked 
the conservation of the Vlasov equation invariants as the norm, the entropy 
and the total energy for the drift-kinetic model. Even if these quantities permit 
to compare the behaviour of the schemes, they do not provide a quantitative 
way to tell one scheme is better than another. We propose an attempt for a 
quality factor Q , based on the dissipation measured with the L'^ norm and 
spurious oscillations measured with the total variation TV norm. It gives a clear 
answer the question which is the best of two schemes, but obviously it is still 
subjective because of the choice of quality factor itself. However, comparison of 
the schemes results with this quality factor is in good agreement with intuition 
when looking at results graphs in ID, but in 4D it helps to evaluate the balance 
between diffusion and spurious oscillations. At the end of this work, we will 
not conclude that one scheme is better than another, we just have given these 
results and tools to help the reader to decide by himself. 
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7 Annexes 

7.1 Annexe A: Detailed LAG reconstruction 

We give hereafter detailed calculations for the LAG reconstruction. 
First, we build the primitive function G at the characteristic feet: 



G{x*j^-yl^ is given by the interpolation with Lagrange polynoms: 

i+2 

G{x) = X] Gj_i/2ij(a;),a; e [a;i_i/2,a;i+i/2] 

j=i-l 

where are the Lagrange polynoms defines as: 

i+2 

^jW = — 



Hence, 



i+2 
k=i—l,k^i—l 



k=i—l,k^i 
i+2 

Li_,{x) = n 

k=i—l,k^i-\-l 

i+2 

L^^ii:^) = n 

fe=i-l,fe#i+2 ■ 

We note x = Xi^i/2 — a,a> then 



a;i-3/2 


- Xk-1/2 


x — 


Xk-1/2 


Xi-1/2 


- Xk-1/2 


X — 


Xk-1/2 


Xi+1/2 


- Xk-1/2 


X — 


Xk-1/2 



i+2 

Li-i{x)= n 

k—i—l,k^i—l 



k^i—l^k^i 
i+2 

.-i(^) = n 

k=i-l,k^i+l 
i+2 
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fe=i-l,;c5^i+2 



2^i-3/2 ^ Xj^_ 
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Assuming Aa; = a;j+i/2 - a;j_i/2,Vi, 



r , ^ TT -a + Axfi - A: + 1) 

= n Ax(^ - 1 - fc) 

r , ^ TT -a + Aa;(z - A: + 1) 

n — A^(— — 

r , TT -a + Aa;(i - + 1) 

= n Ax(i+i-fc) 

-a + Ax(z-fc + l) 
dU.A^)= n A.(i + 2-fe) 



-^+(i-A; + l) 



Writing /? = a/ Ax, we deduce: 

r I \ TT -A^ -r V' - -r 

n --^^^ 

J ( \ TT -^ + (i-fc + l) 

= n (^+2-fc) 

and then, 

r . X _ -/3 + - (i) + 1) + - + 1) + 1) -/3 + - + 2) + 1) 
Li{x) = 



Li+i{x) 

Li+2{x) 



(i-l 


-ii)) 
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-(^ + 2)) 
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- (i + 2) + l) 




1)) 






('; + 2)) 


-f3 + {i- 




-/3 + - (^ + 2) + 




- - (0 + 1) 




-(»-!)) 


(i + l-(i + 2)) 








-(i-l) + l) 


-/3 + (i-(i) + l)- 




-(i + i) + i) 



(i + 2 - (i - 1)) {i + 2- (i)) {i + 2-{i + 1)) 

-/3 + 1 -/3 -/3 - 1 



Li_i(a;) 



1 -2 -3 

L,(x) = ^— — — ^ 

/3 + 2 -;3 - 1) -/3 + 1 



2 -1 1 

+ 2 -/3 + 1 
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L,_i(a;)= 1/6(^-1) (^)(/3 + l) 

L,(a;) = l/2(/3-2)(/3)(/3 + l) 
Li+i(x) = l/2(/3-2)(/3 + l)(,/3-l) 
Li+2(x) = -l/6(/3-2)(/3-l)(/3) 



Li_i(x) = l/6(/33-/3) 
L,(x) = -l/2(/33-/32_2/3) 
Li+i(x) = l/2(/33_ 2/32-/3 + 2) 
Li+2(x) = -l/6(/3-''-3/32 + 2;3) 

Hence 

G{x) = l/6G,_3/2(/3' - /3) 

- l/2G,_i/2(/33 - /32 - 2/3) 

+ l/2G,;+i/2(/^-^- 2/3^-/^ + 2) 

- l/6G,+3/2(/3'- 3/3^ + 2/3) 

Finally, the interpolation can be written: 

G{x) = /33(l/6G,_3/2-l/2G,_i/2 + l/2G,+i/2-l/6Gi+3/2) 
+ (l/2Gi_i/2 - G,+i/2 + l/2G,+3/2) 

+ ^(-l/6Gi_3/2 + Gi_i/2 - l/2Gi+i/2 - l/3Gi+3/2) 
+ ^^1+1/2 

We note x = Xi^i/2 — Aa;/3 = Xi_i/2 — Ax{0 — 1) 
and by defining 6 = 1 — P>0, x = Xi_i/2 + Ax6, then 

G{x) = l/6G,_3/2((l - - (1 - ^)) 

- l/2G,_i/2((l - ef - (1 - ef - 2(1 - 6)) 

+ l/2Gi+i/2((l - Of - 2(1 - ef -{1-0) + 2) 

- l/6G,+3/2((l - Of - 3(1 - 0f + 2(1 - 6))) 

G{x) = l/6Gi_3/2{-0^ + ^0^-29) 

- l/2Gi_i/2i'0'' + 29^ + e-2) 
+ l/2G,+,/2{-e^ + 6^+26) 

- l/6Gi+s/2{-0' + e) 

Gix) = ^3(-l/6G,_3/2 + l/2Gi_i/2 - l/2Gi+i/2 + l/6G,+3/2) 
+ ^^(l/2Gj_3/2 - Gi_i/2 + l/2Gj+i/2) 
+ 9{-l/3Gi_3/2 - l/2Gi_i/2 + Gi+1/2 - l/6Gi+3/2) 

+ Gs_i/2 
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Using the previous formula with x = a;j_|_i/2 — Aa;/3, /3 < 0, i ^ i+l, — >■ — /3, 
we obtain for negative displacement 

G{x) = /33(l/6Gi_i/2 - l/2Gi+i/2 + l/2Gi+3/2 - 1/6^^+5/2) 

+ /?2(l/2G,_i/2 - G,+i/2 + l/2G,+3/2) 

+ /?(l/3G,_i/2 + l/2Gi+i/2 - Gi+3/2 + l/6Gi+5/2) 

+ Gi_|_i/2 

We obtain the following results with x = x^^i/^ — Ax/3: 
if ^ > 0, 

G{x) = /33(l/6G,_3/2-l/2G,_i/2 + l/2G,+i/2-l/6Gi+3/2) 

+ /32(l/2G,;_i/2 - G,+i/2 + l/2G,+3/2) 

+ /3(-l/6G^_3/2 + Gi_i/2 - l/2Gj+i/2 - l/3Gj+3/2) 



if /3 < 0, 

G{x) = /33(l/6Gi_i/2 - l/2Gi+i/2 + l/2Gi+3/2 - l/6Gi+5/2) 

+ /32(l/2Gi_i/2 - Gi+1/2 

+ l/2G,+3/2) 

+ /3(l/3G,_i/2 + l/2Gi+i/2 - Gi+3/2 + l/6Gi+5/2) 

+ Gj+1/2 

We define gfAx = Gj+1/2 — Gi_i/2- The previous expressions become: 

if /3 > 0, 

G{x) = /33(l/6Gi_3/2 - l/6G,_i/2 - 2/6G,_i/2 

+ 2/6G,+i/2 + l/6G,+i/2 - l/6G,+3/2) 

+ /32(l/2Gi_i/2 - l/2G,+i/2 - l/2G,+i/2 + l/2Gi+3/2) 

+ /3(-l/6Gi_3/2 + l/6G,_i/2 + 5/6Gi_i/2 

- 5/6Gj+i/2 + l/3Gj+i/2 - l/3Gj+3/2) 

+ ^^1+1/2 

= /33(-l/65r-i + l/S.gf - l/65?+i) 

+ /32(-l/25r + l/2.9r+i) 

+ /3(l/65r-i - 5/6<?r - l/3<??+i) 
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if /3 < 0, 

G{x) = /?3(l/6G,_i/2-l/6G,+i/2-2/6G,+i/2 

+ 2/6Gi+3/2 + l/6G,+3/2 - l/6G,+5/2) 

+ /3'(l/2G,_i/2 - l/2G,+i/2 - l/2G,+i/2 + l/2G,+3/2) 
+ /?(l/3G,_i/2 - l/3G,+i/2 + 5/6G,+i/2 

- 5/6Gi+3/2 - l/6Gi+3/2 + l/6Gi+5/2) 
+ G i+1/2 

= /33(-l/6<7r + 2/65r+i-l/65?+2) 

+ /32(-i/25r+i/25r+i) 

+ /?(-i/3gr-5/65r+i + i/6flr+2) 

+ Gi+i/2 

Let now compare theses fluxe expressions with the regular expression we 
found for the Hermite formulation noticed (t>k+i/2,jk+i/2=k{l^)- 

7.2 Positive displacement: /3 > 0,5 — l,jk+i/2 — k 

Sk (^'(3-2/3)) 
/3^(.9fe-i/2+fffe+i/2-25fe) 
/^'(-.9^-i/2 -2fffe+i/2 + 35fc) 

/3(5fe+l/2) 

We have introduced the value of the distribution at the cell's faces that we note 
9k-i/2 = dix ^ Xk-1/2) and 5^+1/2 = 9{x ^ Xk+i/2)- We set these coefficients: 

j9ti/2 = W-i + 5/65^ - 1/65^+1 
l/6ffLi + 5/65^ + l«+i 



9k~i/2 
<7fc+i/2 

5fe (/3'(3-2/?)) 

^3(l/35^_i + 5/65;' - I/69U1 - 1/65^-1 
5/6ff^ + l«+i-25^) 
/32(-l/35^_^- 5/65^ + 1/65^+1 
l/3.9Li - 10/6.9^ - 2/3.9,"+! + 3<?^+i) 
/3(-l/65fe"_i + 5/6<?^ + l«+i) 



'?^fc+l/2,ifc+i/2=fc(/3) 

Aa; 



+ 
+ 



+ 



{9k+l/2 

that leads to: 

0fc+l/2,jfc + i/2=fc(^) 

Aa; 



+ 
+ 
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- .9fe-l/2 [P [P - l)j 

+ fffe+1/2 - /3f ) 

+ 5fe(/3'(3-2/3)) 

= ^3(l/65Li-l« + l«+i) 

+ P\l/29l-l/29l^,) 

+ ^(-l«_i+5« + l«+i) 

Remembering that 

= ^3(-l/65?_i + - l/65?+i) 

+ /32(-l/2gr + l/25r+i) 

+ /3(l/6ffr_i - 5/65? - l/35?+i) 

we deduce: 

G{x) = -0fc+i/2jfc+i/2=fe(/3) + G,+i/2 
7.3 Negative displacement 5 = 

'/'fe+l/2,Jfc + i/2=fe + l(/3) 



Aa; 



= 5fc+l/2 (W + 1)') 

+ 5fc+3/2(/3'(l + /3)) 

+ 5^+i(/32(-3-2/?)) 

= /3^(5fc+3/2 + fffc+l/2 ~ 25^+1) 

+ /3'(25^+i/2+flfc+3/2-35fc+i) 

+ /3(5^+l/2) 



G(x) = /33(-l/6ff? + 2/6<?r+i-l/65?+2) 

+ /32(-l/25? + 

+ /3(-l/3ffr - 5/6ffr+i + l/65?+2) 

+ Gj+1/2 



We set these coefficients: 



\9t+i/2 = ^IM + 5/65^+1 - 1/65^+2 

[9k+z/2 = -1/65^ + + 1/35^+2 



0fc+l/2,jfc+i/2=fc+l(/3) 



Ax 



+ 5/6ff^+i + l/3ff^+2 - 25^+1) 

+ /32(2/35^ + 10/65fc"+i - 2/6g^+2 - 1/65^ 

+ 5/65^+1 + 1/35^+2 - igl+i) 

+ /3(l« + 5/6<?^+i- 1/65^+2) 
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Aa 



/?^(l/6.9l?-l/35lVi + l/6ff^2) 
+ /?(l/3fffc" + 5/6ff,V-l/65;V2) 



Remembering that: 



G{x) = /33(-i/65r + l/3.gr+i-l/6gr+2) 
+ /32(-l/25r + l/2gr+i) 
+ ^(-l/a^r - 5/65r+i + l/6gr+2) 

we deduce 

G{x) = -<?ifc+i/2j,^^/,=fc+i(/3)/Ax + G,+i/2 

To conclude, 

ffr" - -'^fc+i/2j.+,/. (/3)/Ax + G„+i/2 + Gfc_i/2^,,_,,, (/3)/Ax - G,_i/2 
ffr"' = 97 - (-^^+1/2,,.,,/. (/?) + <Pk-i/2,j._,,. (/3))/Ax 

7.4 Annexe B: Algorithms 

We present some algorithms which are representative of the Gysela algorithms. 
We change the notation of the flux 0i+i/2 to i?i+i/2 for practical reasons. The 
ID advection equation are solved with Algo. ((1]). 
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Algorithm 1 Advection 

Require: (a;fc)fc=-i/2...w+i/2 > i9k)k=o-N 
1; a=l,b=-l,c=0,p=l,q=4 
2: L,U=MatrixCoefRcient(a,b,c,p,q) 

3' {9k-i/2)^_^ ^^^==CoefRcientHermite((5^)^^o...^,L,N,0,N+l,a,b,c) 

4: /3=FootCharacteristic(...) 
5; if /3 = then 

6: -ffo-1/2 = 

7: else 

8; if ^ > then 
10: else 

11: ffo-i/2=FluxComputation(gJ^, 5^_^/2' 5'o+i/2> ^) 

12: end if 

13: end if 

14: for i=l:N-l do 

15: /3=FootCharacteristic() 

16: if /? = then 

17: Hi+i/2 = 

18: else 

19: if /3 > then 

20: j=i 

21: else 

22: j=i+l 

23: end if 

24: i?j+i/2=FluxConiputation(5j', ff"_i/2. 5"+i/2. /3) 
25: end if 
26: end for 

27: j5=FootCharacteristic(...) 
28: if = then 

29: Hn+1/2 = 

30: else 

31: if /? < then 

32: Hn+1/2 = PQn 

33: else 

34: i?Ar+i/2=FluxComputatioil(ff]^, fif]^_i/2' 5]^+l/2' P) 

35: end if 
36: end if 

37: for i=l:N-l do 

38: 5?+' = - {H,+l/2 - if(.-l) + l/2) 

39: end for 

40: return Sffc=o-JV 
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7.4.1 FluxComputation and FootCharacteristic 



Algorithm 2 FluxComputation 
Require: g^, 9l_y^, g'l+y^, f3 
1; if /3 > then 

2: 5 = 

3: else 

4: 6 = 1 

5: end if 

6: i?fc+i/2 =3,.-i/2 (/3(l-<5)+/32(2-3(5)+/33) 

7; = i?fc+l/2 + 5,. + l/2 + - 35) + 

8; Hk+1/2 = Hk+1/2 +9j, (/32(-3 + 65)+/33(-2)) 
9: return i?fc+i/2 



Algorithm 3 FootCharacteristic 
1: return x* , ,„ 

1+1/2 



7.4.2 Hermite coefficient computation for natural conditions 

For a domain with natural boundary conditions, we compute the value of the 
distribution at the cell's faces gk-1/2 and gk+i/2- We want to solve: AX = B 
and we decompose A &s A = LU with Algo. ([5]) then we solve X = U^^L^^B 
with Algo. (gl). 

Algorithm 4 CoefficientHermiteNatural 

Require: {gD^^^.^j^ , {Lk)k=o-N+i ' {Uk)k^o-N+ifi^^ + l'a,b,c 
1: Xo = 5gE 

for fc = 1 : iV do 

Xk = 3{g^_,+gl)-Lk-iXk-i 
end for 



7: 
8: 
9: 
10 
11 
12 



Xn+1 — 5.g]^ — Ln+iXn^i — L^X^ 

V ^N + l 

^N+l - 

for k = N -.2, 1 do 

Xk 
end for 

_ Xi-{l-c/a)X2 



^fe - — m — 



return (.gfc-i/2);,=o. 



y _ Xo-bXi-cX2 
^0 - Uo 



■N+1 



We compute the matrix coefficients of the LU decomposition. 

[A] = [L] [U] 
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Algorithm 5 MatrixCoefficientNatural 



Require: a,b,c,d,e,f,p,q,0,N+l 



1 


Uo = a 




2 


^0 = ^ 




3 


Ui = q- 




4 


Li = i: 




5 


c/2 = 4- 




6 


for k — 


2 : TV- 1 


7 


Lk = 




8 


Uk+i 


= q - Lk 


9 


end for 




10 


Ln+1 = 


d 


11 


Ln — 


-^N + l 

Un 


12 


dN+i = 


f-UN 


13 


return 


L,U 



7.4.3 Hermite coefficient computation for periodic conditions 

We compute the value of the distribution at the cell's faces gk-1/2 ^nd gk+i/2- 
With periodic boundary conditions, we decompose A &s A ~ LDL^ with Algo. 
then we solve Y = D~^L~^B with Algo. 
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Algorithm 6 CoefRcientHermitePeriodic 



Require: ((?^)fc^o...jv_i , {Lk)k=o-N > (^^^fe)fe=o...iV'a,b,c 
1: Xo = 3(5]^ + 5o") 
2: for = 1 : AT - 1 do 
3: Xk = 3{g^ + g^^,) -Lk-iXk-i 
4; end for 

5: S=0 

6: for m = : iV - 2 do 

7; 5 = 5* + 6mXm 

end for 

Xn = + .9o ) - sum - Ln-iXn-i 
for m = : N do 



9 
10 
11 
12 
13 
14 
15 
16 
17; 



Y — 2£nL 
^™ ~ Dm 

end for 

Xn-i ~ Xjy^i — L(^N — 1)Xn 
for fc = AT - 2 : 0, -1 do 

Xk = X}: — LhXk+i — SkXjsr 
end for 

return {gk-1/2) ^^o-n 



We compute the matrix coefficients of tlie LDL decomposition. 

[A] = [L] [D] 
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Algorithm 7 MatrixCoefRcientPeriodic 



Require: start -.i end 
Di = q 

for fc = 2 : iV - 2 do 
Dk=q- pLk-i 

T, P_ 

end for 

Dn-1 = q-pLN-2 

J- p-pSn-2 

DN = q- Ef=i' Di6,^ - Dn-iL%_ 
return L, D 
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